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The Kolmogorov-Johnson-Mehl-Avrami (KJMA) theory 
for the time evolution of the order parameter in systems un- 
dergoing first-order phase transformations has been extended 
by Sekimoto to the level of two-point correlation functions. 
Here, this extended KJMA theory is applied to a kinetic Ising 
lattice-gas model, in which the elementary kinetic processes 
act on microscopic length and time scales. The theoretical 
framework is used to analyze data from extensive Monte Carlo 
simulations. The theory is inherently a mesoscopic continuum 
picture, and in principle it requires a large separation between 
the microscopic scales and the mesoscopic scales characteristic 
of the evolving two-phase structure. Nevertheless, we find ex- 
cellent quantitative agreement with the simulations in a large 
parameter regime, extending remarkably far towards strong 
fields (large supersaturations) and correspondingly small nu- 
cleation barriers. The original KJMA theory permits direct 
measurement of the order parameter in the metastable phase, 
and using the extension to correlation functions one can also 
perform separate measurements of the nucleation rate and 
the average velocity of the convoluted interface between the 
metastable and stable phase regions. The values obtained for 
all three quantities are verified by other theoretical and com- 
putational methods. As these quantities are often difficult to 
measure directly during a process of phase transformation, 
data analysis using the extended KJMA theory may provide 
a useful experimental alternative. 

PACS numbers(s): 64.60.Qb, 75.60.-d, 77.80.-e, 82.65.-i 



I. INTRODUCTION 

The phase-transformation kinetics in systems undergo- 
ing first-order phase transitions are important in many 
scientific and technological disciplines. Around 1940, 
Kolmogorov,El Johnson and Mehl,H and AvramH (KJMA) 
introduced a simple theory describing the decay of a 



metastable system towards a unique equilibrium phase. 
This theory applies to systems with a nonconserved order 
parameter, in which the decay is driven by a difference 
between the free-energy densities of the metastable and 
equilibrium phases. Originally conceived for metallurgi- 
cal applications, it was later formalized and generalized 
by Evans,cl who also pointed out its applicability in sur- 
face science. 

The basic assumptions of the KJMA theory are sim- 
ple: negligibly small "droplets" of the equilibrium phase 
nucleate from the uniform metastable phase and sub- 
sequently grow without substantial deformation. The 
growing droplets are assumed to be randomly placed and 
overlap freely, with the result that the remaining volume 
fraction occupied by the metastable phase decays expo- 
nentially with a power of time: 



<p ms (t) =exp(-Ct tt ) . 



(1.1) 



Both the coefficient C and the "Avrami exponent" a de- 
pend on the spatial dimension d and on details of the 
nucl eatio n and growth processes. In general, a > 1; thus 
Eq. (1.1) does not represent a stretched exponential. 

The KJMA picture has been extensively applied 
in diverse fields of research.0 A small selection of 
examples include transitions between different liquick 
crystal.-. phases;tl .crystallization kinetics p._Lipids,Ll 
sugarsja polymers,Ei and eute^ic-|mixtures;E30 solid- 
state phase tpj^>£prmations;lljilj dotmaia switching 
in ferroelectricsEirEj and ferromagnets;Ejc3 adsorption 
and surface growth kinetics, in electrochemicallffij 42 ! or 
gas/vacuum environments pj rock format ion ;E3cJ and 
slow combustionJla 

During the first forty years after its inception, the 
KJMA theory appears to have been appreciated mostly 
by experimentalists. It attracted little sustained theoret- 
ical interest until the 1980s, when Sekimoto derived ex- 
act expressions for the two-point-phase correlation func- 
tion within the KJMA pictureJl3 Knowledge of corre- 



1 



a fixed shape nucleate, either deterministicall 
domly, and thereafter grow deterministically. 



lation functions makes it possible to predict and inter- 
pret the results of small-angle scattering experiments. In 
this paper we demonstrate how it can also be utilized to 
obtain independent estimates of the droplet nucleation 
rate and growth velocity. Sekimoto's results have been 
generalized to .-systems with infinitely degenerate equi* 
librium phasesli3 and to the case of finite degeneracy.ED 
We shall refer collectively to the theories that extend 
the KJMA picture to include correlation functions as 
the "extended KJMA theory." Among recent applica- 
tions of the extended KJMA theory arc .theoretical stud- 
ies of domain switching in ferroelectricsta~c3 predictions 
of magnetic-fpKC|e-JTiicroscopy observables for nanoscale 
ferromagnetsJ^SnS^LsLftd^tudies of hysteresis in spatially 
extended sy stems. t3Ejc§c2l Relations between two-point 
correlation functions and droplet size distributions have 
also been discussed recentlyEj 

There are many simplifying assumptions inherent in 
the extended KJMA theory. These include a constant nu- 
cleation rate, a constant interfacial propagation velocity, 
and the neglect of surface-tension effects that one would 
expect to become important when growing droplets meet 
and coalesce. In view of this, it is remarkable that the 
theory appears to perform as well as it does for such 
a wide variety of systems. Computer studies have pre- 
viously been performed for models in which droplets of 
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However, we are only aware of a single previous studyu in 
which the kinetic processes act on length and time scales 
that are microscopic compared to the growing droplets, 
as one would expect for real physical and chemical sys- 
tems. Such studies are desirable, as they permit direct 
verification of the regimes of validity for the individual 
assumptions as well as the consequences of particular as- 
sumptions not being exactly fulfilled. In the present pa- 
per we further address this need, with particular empha- 
sis on the information which can be obtained by com- 
bining results for the volume fraction and the two-point 
correlation functions. Our purpose is two-fold, as dis- 
cussed below. 

Our first aim is to perform a detailed test of the valid- 
ity of the extended KJMA theory in a particular model 
system with microscopic kinetics. We consider a two- 
dimensional kinetic Ising model of ferromagnetic or ferro- 
electric switching. This is equivalent to a simple lattice- 
gas model of submonolayer chemisorption onto a single- 
crystal surface in the limit that lateral diffusion can be 
ignored.c3 For this model we test the KJMA predictions 
for time-dependent quantities that describe the meso- 
scopic two-phase structure during the evolution towards 
equilibrium. These quantities include the volume frac- 
tion, the two-point correlation function, and its Fourier 
transform, the structure factor. We find excellent quanti- 
tative agreement between the theoretical predictions and 
the simulation results for a considerable range of applied 
fields and times. 

Our second aim springs from the significant regime 



of agreement that we establish between the KJMA pre- 
dictions and the simulation results for the model stud- 
ied. This enables us to use the KJMA prediction for 
the volume fraction to subtract those regions of the sys- 
tem which have already decayed to the stable phase 
at a particular time, and thus measure extended time 
and volume averages of the quasi-equilibrium order pa- 
rameter in the metastable phase. This approach rep- 
resents a novel method to measure "thermodynamic" 
quantities in a constrained metastable ensembleO Fur- 
thermore, the extended KJMA predictions for the two- 
point correlation functions enable us to measure sepa- 
rately the nucleation rate and the average propagation 
velocity of the convoluted, moving interface between the 
metastable and stable phase regions. We show that the 
measured quantities agree well with theoretical predic- 
tions obtained by different methods. These methods are 
numerical transfer^ matrix calculations for the metastable 
order parametecXTO a field-theoretical result for the nu- 
cleation rate,E2lEil and— a nonlinear response theory for 
the interface velocityj£3 The close agreement between 
the KJMA and independent estimates confirms that the 
KJMA theory does not merely produce good fits to the 
simulation data, but that the fitting parameters measure 
nonequilibrium physical quantities that can be difficult 
to measure by other methods. 

The remainder of this paper is organized as follows. 
In Sec. H we introduce the model and the numerical 



methods used in this work. In Sec. [II 



we summarize 



those results of the nucleation theory of metastable decay 
and the extended KJMA theory that are relevant to our 
study. We give most of the theoretical results for general 
spatial dimension d. In Sec. IV our two-dimensional sim- 



ulation results are presented and explained in the frame- 
work of the theory. The regimes of agreement between 
theory and simulations are identified, and the simulation 
data are used to measure the metastable order parame- 
ter (Sec. [VA), as well as t he ave rage interface velocity 
and the nucleation rate (Sec. IV B ), all a s func tions of the 
applied field. Co rrelat ion functions (Sec. IV C) and struc- 
ture factors (Sec. [TV D| ) are also measured. The measured 
quantities are compared with independent theoretical es- 
timates. In Sec. |v| we summarize our results and give our 
conclusions and some suggestions for further study. 



II. MODEL AND METHODS 

A. Model Hamiltonian 

The discrete model studied here is defined by the stan- 
dard Ising Hamiltonian (the transformation to lattice-gas 
language is given below), 



H = -J^Si(t) Sj(t) -H^2si(t) , 



(2.1) 



(i,3) 
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where Si(t) = ±1 is the "spin" at lattice site fi at time 
t, J > is the ferromagnetic coupling constant, and 
H is the external field. The sums Yltij) an d Si run 
over all nearest-neighbor pairs and over all sites on a 
rf-dimensional hypercubic lattice of linear size L, respec- 
tively. The lattice constant defines our unit of length. 
While our theoretical discussion is for general d, all the 
simulations presented are for d=2. Three-dimensional 
systems will be considered in a forthcoming paper S In 
order to avoid complications due to boundaries, we use 
periodic boundary conditions throughout. (For recent 
discussions of boundary effects in metastable decay, see 
Refs. |9]|4].) 

The magnetization per unit cell is 



(2.2) 



Under equilibrium conditions m is the order parameter 
conjugate to H. For H=0 and temperatures below a fi- 
nite critical temperature T c , the model has two degener- 
ate equilibrium phases in which the magnetization has a 
constant spontaneous magnitude m sp (T) = m s (T, H—0). 
In the presence of a nonzero field the phase degener- 
acy is lifted, and the stable equilibrium magnetization 
m s (T,H) has the same sign as H. 

Although it is not absolutely essential for our study 
that T c and m sp (T) are exactly known ..for the two- 
dimensional square-lattice Is ing m odel, E3i£l3 these and 
other exact results (see Sec. Ill A) enable us to quan- 
titatively compare the various parameter estimates ob- 
tained from our numerical simulations with independent 
theoretical estimates. 

The Ising formulation is conveniently symmetric un- 
der simultaneous reversal of H and Si(t), and it can be 
directly applied as a simple model for highly anisotropic 
ferromagnetic and ferroelectric systems. A less explic- 
itly symmetric formulation is particularly convenient 
for discussing crystallization and adsorption phenom- 
ena. Itis-the equivalent two-state, attractive lattice-gas 
model.OH In terms of the time-dependent, local occu- 
pation variables c;(i) € {0, 1}, Eq. (2.1) takes the form 



L d ( 1 



(2.3) 



The quantities appearing i n th e equi vale nt formulations 
of the Hamiltonian, Eqs. (2.1) and (2.3), are linked by 
the transformations 



Ci(t) = \ [1 



4 J , 
2H 



Mo 



(2.4a) 

(2.4b) 
(2.4c) 



Here $ is the attractive lattice-gas interaction energy, 
and \i is the chemical or electrochemical potential, whose 



value at coexistence (i.e., for H=0) is fio = —2zJ = 
— z<&/2, where z is the coordination number (z = 2d for 
hypercubic lattices). The chemical potential is related to 
the (osmotic) pressure p as fi — fj,Q = fcsrin(p/po), where 
Pq is the pressure at coexistence, T is the temperature, 
and fee is Boltzmann's constant. The order parameter 
conjugate to [i is the density (for d=2: the coverage), 



9(t) 



(2.5) 



B. Stochastic Dynamics 

The Ising lattice-gas Hamiltonian does not impose a 
particular dynamic on the system. To study the ap- 
proach to equilibrium under the influence of thermal fluc- 
tuations, we use the stochastic Glauber dynamicO This 
dynamic is defined by the acceptance probability for a 
proposed flip of Sj, 



W\ 



exp(-/3AE) 
1 + exp(-/3A£) 



(2.6) 



where AE is the energy change associated with the at- 
tempted spin flip, and (3 — l/k^T. The dynamic was 
implemented by the standard algorithm used for kinetic 
studies:L°J a series of steps in which sites are chos en a t 
random and flipped with probability given by Eq. (2.6). 

The Glauber-Ising model described above (or the 
same Hamibxinian with the closely related Metropo- 
lis dynamicLil) has previously been used to study a 
large number of phase-ordering phenomena in condensed- 
matter physics and other fields. The two-dimensional 
version corresponding to the model for which we 
present numerical results, has been applied to Acfc 
mally activated sadtching in uniaxial ferroelectricE-J~t3 
and ferromagneticE3~E3 thin films. The equivalent lattice- 
gas model should give a reasonable representation, pf 
the kinetics of-submonolayer adsorption in chemicalEj c3 
and physicalcj systems in which lateral adsorbate dif- 
fusion canJae ignored. Electrochemical underpotential 
deposition^ in which second-layer formation is heavily 
suppressed, might be approximated by this model,E3 al- 
though more detailed agreement with experimental stud- 
ies of dynamical phenpmcnaEil requires the inclusion of 
lateral diffusion.BEO 



C. Quantities Characterizing the Decay 

We study the decay of a metastable phase by first 
preparing the system in an initial state of magnetiza- 
tion m(0)—+l. We then apply a constant magnetic field 
H < and let the system evolve in time. The decay 
of the metastable phase is characterized by the behavior 
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of the relaxation function,0 which is defined in terms of 
m(t) as 



m(t) — m s 
m(0) — m s ' 



(2.7) 



where m s (T, H) is the field-dependent equilibrium mag- 
netization. For m s = — 1, <j){t) simply e qua ls the time- 
dependent density 9(t) defined in Eq. (2.5). The life- 
time t of the metastable phase is defined as the average 
first-passage time to to=0. In the parameter regime de- 
scribed by the KJMA theory, this is equivalent to the 
requirement that the ensemble average (m(r))=0. (For 
a discussion of the effects of using different cutoff values 
of m to define r, see Ref. [75} ) 

During the decay process, we also compute the cir- 
cularly averaged time-dependent structure factor S(q, t) 
and correlation function G(r,t). The time-dependent 
structure factor S(q, t) is defined in terms of the Fourier 
transform of Si(t), 



-iq-ri 



(2.8) 



as 



S(q, t) = (%(t) ~ M*))\o ■ ( 2 - 9 ) 

The components of the reciprocal lattice vectors q are 
q ja = 2wj/L [j = 0,±l,±2,...,±(L/2- l),L/2; a = 
x,y, ■ ■ .], and is the Kronecker delta function. The 
brackets ( ) imply an ensemble average over independent 
simulation runs. 

The time-dependent two-point correlation function 
G(r , t) is defined by 

G{r,t) = ( Si (t) 8 S (t)) - (m(t)) 2 , (2.10) 

where r=rj — fj. It is circularly symmetric an d wa s ob- 
tained as the inverse Fourier transform of Eq. (2.9). All 
Fourier transforms were computed with the Fast Fourier 
Transform subroutine fourn.L3 With the normalizations 
used here, S(q = 0, t) = L d Var[m(i)], which is indepen- 
dent of L if G(r,t) is of finite range. Here, Var[m(t)] 
is the time-dependent variance of the system magneti- 
zation, evaluated over an ensemble of independent runs. 
The quantities S(q, t) and G(r, t) were circularly aver- 
aged to obtain S(q,t) and G(r, t), respectively. 

In order to find a time-dependent characteristic length 
scale, the first moment of G(r, t) was computed using the 
definition 

LrG(r,t) 



(K*)> 



J2 r G(r,t) 



(2.11) 



The temperature in all our simulations was fixed at 
T = 0.8T C , which is sufficiently far below T c that the ther- 
mal correlation lengths in both the stable and metastable 
phases are microscopic and the main contributions to the 
correlation function come from the random two-phase 
structure. Other details about the simulations are given 
at the beginning of Sec. IV. 



III. THEORETICAL BACKGROUND 
A. Nucleation Theory of Metastable Decay 

In order to obtain the characteristic times and lengths 
that are important to our analysis, we here give a brief 
summary of those aspects of the droplet theory of nucle- 
ation that are necessary to analyze our numerical results. 
More complete discussions can be found in Refs. p4] , [75| , [77| . 

Thermal fluctuations in the uniform metastable phase 
create droplets of the stable phase. In terms of the 
droplet radius R, the free-energy difference between a 
uniform metastable system and one that contains a sin- 
gle such droplet is 

AF(R) = dnp^R d - 1 cTo - fl d R d \H\ ■ \m s - m m8 | , (3.1) 

where QdR d is the volume of a droplet of radius R, <tq 
is the surface tension along a primitive lattice vector, 
and m ms is the magnetization of the uniform metastable 
phase. The critical radius 



R c = 



(d - l)o-o 



(d-l)o-o 



H 



2\H\m 



(3.2) 



sp 



corresponds to the maximum of AF(R). Droplets with 
R < R c almost always decay, whereas droplets with 
R> R c almost always continue to grow. In this work we 
only consider sufficiently strong fields that R c <SiL. (The 
regimes of very weak fields, where this does not hold, are 
discussed in Refs. |^-|9|,|4|,^5U7|.) The nucleation rate 
per unit volume for growing regions of the equilibrium 
phase is determined by the free-energy cost of a crit- 
ical droplet, AF(R C ), through a Van't Hoff-Arrhcnius 
relation. J^ bc ,. fi nal form is shown by field-theoretical 
argumenta_3'EaEil to be 



I(T,\H\)nB(T)\H\ K exj>(- 



E(T) + 0(H 2 ) 



\H 



Id— 1 



(3.3a) 



where 



E(T) 



1 



k B T Vl?™ s -ra ms | 

n d a d U-V "' 1 
k B T \2m sp 



d-l 



(3.3b) 



r. There 
9 ' evidence 



In Eq. (3.3a), B(TX.is-.& nonuniversa' 
is strong analyticalE£lEll and numericaSE 
that the exponent if is 3 for the two-dimensional Ising 
model, and it is believed-, to be —1/3 for the three- 
dimensional Is ing m odel. E3 The approximate form of 
S(T) in Eq. (3.3b) is completely defined by quanti- 



ties that for the two-dimensional f.sing model are either 
known exactly (cq and m sp ),EJEj or can be obtained 
through a Wulff construction by pnmeximl integration 
of exactly known quantities (n d )BBEEl The Q(H 2 ) 
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corrections in Eq. (3.3a) are relatively minorB and will 
be ignored. 

The above calculations are based on the continuum 



assumption that 2R C ^> 1. From Eq. (3.2) one sees that 
this requires that 



H « (d - l)a (T)/m sp (T) = F M fsp(T) 



(3.4) 



This crossover field has been called the mean-field spin- 
odal (MFSP), and the regime of stronger fields is called 
the strong-field regimeOH Note that -Hmfsp depends 
on T, but not on L. 



B. Continuum KJMA Theory 

1. One-point averages 

The KJMA theory of metastable decay describes the 
process of nucleation and growth in a large continuum 
system with a nonconserved order parameter and a non- 
degenerate equilibrium phase. (The precise meaning of 
"large" will be elucidated below.) It is assumed that 
the system is initially in a uniform metastable phase, in 
which critical droplets of the stable phase nucleate with 
rate I(t) per unit volume and grow with radial growth 
velocity v(t) without substantially altering their shapes. 
The two-phase structure of the system is represented by 
a phase field or indicator function, 

« J 1 if r is in the metastable phase at t , . 
^ ' ' [0 otherwise . \ ■ ) 

The simplest statistical quantity describing the structure 
of the system is the volume fraction of metastable phase, 
<Pms(t), which is defined in terms of u(r,t) by 

<Pw,(t) = (u(t)) = -1 f (u(r, t)) d d r . (3.6) 



Assuming that droplets of stable phase nucleate wit h con - 
stant nucleation rate I (the case discussed in Sec. [II A ) 
and grow from an initial volume of zero (R c <C L) with- 
out interacting and with constant radial growth velocity 
v(t) = V, the ensemble average of ip ms (t) is given by 



<Pms(i) = exp 



exp 



-m d v d / (f - sfds 
Jo 

n d iv d t d+1 



d+l 



(3.7) 



an expression often referred to as "Avrami's law."EJ~El 
The argument of the exponential (which grows with- 
out bound with time) is the "extended volume" of sta- 
ble phase, obtained by adding the volume fractions of 
all droplets without correcting for overlaps. The expo- 
nential dependence on the extended volume is a special 
case of a general result for randomly placed, polydis- 
perse objects. In the approximation considered here the 



distribution of droplet radii is uniform between R c « 
and vt. In response to claims that it does not repre- 
sent a correct solution for the stochastic process defined 
in Refs. Eq. (3__7) was recently ijederived without 
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reference to the the extended volumea In Sec. 
Eq. (|3.7|) will be used to provide explicit form s of the pa- 
rameters in Eq. (fLl|). Generalizations of Eq. ( |3.7[ ), which 
consider complications such as finite-size effects, homoge- 
neous and heterogeneous nucleation, anisotr opic g rowth 
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and diffusion, are discussed in Refs. |14| , |24[ 
Effects of nonstationary nucleation rates and growth ve- 
locities that depend on the droplet size are discussed in 
Refs |5lM|S7l 



Equation (3/7) defines the time scale to = A (lv d ) d+1 , 

i 

in which A = \{d + 1) In2/Qdj d+1 depends weakly on 
T through fi^. This time approximately equals the 
metastable lifetime r defined after Eq. (2.7). An im- 



portant length scale is obtained from to and the growth 
velocity v: 



Ro = vt = Av (lv d ) ' 



(3.8) 



This characteristic length describes the mesoscopic struc- 
ture of the decaying system. It gives the average diam- 
eter of a droplet at t » r and can be seen as the aver- 
age distance between independent droplets. The average 
number of droplets that contribute to the decay is pro- 
portional to (L/Ro) d . For Eq. (3.7) to describe the time 
evolution correctly, the system must contain a large num- 
ber of independently nucleating and growing droplets, 
i.e., (L/Ro) d ^>l. This is the sense in which the system 
must be large. Because of the large number of droplets 
that contribute to the growth, the regime in which KJMA 
theory is exp ected to be valid is called the multidroplet 
(MD) regime£3uj Under these conditions the system is 
self-averagingE3 and behaves approximately determinis- 
tically according to Eq. (3.7). 

In terms of the four characteristic lengths - the micro- 
scopic lattice constant (unity), the critical droplet radius 
i? c , the average droplet distance Ro, and the system size 
L - the domain of validity of the KJMA approximation 
can be summarized as 



1< R c < Ro < L 



(3.9) 



Setting Ro ~ L, one obtains the crossover field that 
limits the MD regime in the weak-field/smalLgistem di- 
rection, called the dynamic spinodal (DSP)0't3 



i SCO 

+ 1 InL 



(3.10) 



Its exceedingly slow asymptotic convergence with L re- 
sults from Eqs. ( [3.3a ) and ( |3.8| ). However, relatively large 
systems (L JSJ 10 a — 10 4 for d =2) a re required for the con- 
tribution described by Eq. (3.10) to be lar ger than the 
various correction terms (see Fig. 11 of Rcf. [7a). 



5 



This discussion reveals three restrictions that pertain 
to the applicability of the KJMA picture to real systems, 
as well as to the Ising lattice-gas model: 

1. It applies only in the MD regime of large 
systems and/or intermediate fields. In this regime 
the decay proceeds through a large number of 
droplets which nucleate independently at random 
times and positions and subsequently grow to fill 
the system. 

2. It describes nucleation and growth in a coarse- 
grained sense. This means that the results of the 
theory should agree with those of the Ising model at 
length scales much larger than the critical droplet 
radius R c and should disagree at length scales on 
the order of R c and shorter. 

3. It does not take into account interfacial effects, 
except insofar as they determine the nucleation 
rate. When the volume fraction of stable phase 
is large, the dynamics are dominated by droplet 
coalescence, which is accelerated by the interface 
tension. We therefore expect the theory to disagree 
with the simulation results in this late-time regime. 



2. Two-point correlations 

The connected two-point correlation function for the 
metastable phase, T(r, t), was obtained in closed form 
by Sekimoto under the assumption that th&jdroplets are 
(i-dimensional spheres (02 = 7r, =4-7t/3)o 



3. Relations between KJMA and Ising quantities 

Theoretical approximate expressions for the relaxation 
function <j>(t) and the correlation function G(r, t) of the 
Ising model can be derived from the corresponding quan- 
tities in the KJMA theory. The main assumption is that 
for sufficiently late times, when the mean size of the 
droplets of stable phase is much larger than R c , the meso- 
scopic structure of the Ising model resembles the KJMA 
picture. This assumption is also expected to break down 
at late times, when droplet coalescence becomes impor- 
tant. 

The time-dependent magnetization of the Ising model, 
m(t), is approximately given in terms of the volume frac- 
tion in the KJMA theory as 



u(x, t)u(x + r, t)) — (u(x, t)) 2 
(u(t)) 2 {exp [lv d t d+1 ^ d {r/2vt)] 
0, 



-1} 



r < 2vt 
r > 2vt 

(3.11 



where r=\r I and 



Mv) = 3 



arccos y 



-y 3 in 



i + - y 2 



*s(y) = -(l-y) 3 (l + y). 
The first moment of T(r, t) is defined by 
/ r r(r, t) dr 



it)) = 



/ r(r, t) dr 



(3.12a) 
(3.12b) 

(3.13) 



consistent with Eq. ( |2.1l| ). This is a time-dependent 
characteristic length which describes the structure of the 
system. One would expect its value at t = r to be pro- 
portional to Rq, which is confirmed by our simulations. 



i(t) 



Vms(i) 



(3.14) 



where m ms and m s are the magnetizations of the domains 
of metastable and stable phase, respectively. From the 
definition of the relaxation function <fi(t), Eq. (2.7), one 
also has 



m(t) = [m(0) — m s ] <p(t) + m s 



(3.15) 



for the Ising model. Equations ( 3.14 ) and ( |3 .15 ) together 
yield 



m(0) 



■ ¥>ms(*) 



(3.16) 



for the relaxation function cf>(t) of the Ising model in 
terms of (p ms (t) in the KJMA theory. The right-hand 
side of Eq. |Tl§) can be considered as a coarse-graining 
approximation for <fr(t) , in which the local spins have been 
averaged over a region which is large compared to i? c , but 
small compared to Rq. 

At this level of coarse-graining, the correlation function 
of the Ising model is given in terms of the corresponding 
quantity in the KJMA theory, T(r, t), as 



G(r,t) 



T(r,i)+G ms (r)+G s (r) . (3.17) 



Here G ms (r) and G s {r) are correlation functions de- 
scribing local fluctuations that are nonzero only in the 
metastable and stable phase, respectively. (See Ap- 
pendix |a|.) These correlation functions are of very short 
range compared to T(r,t), and where they are different 
from zero they are proportional to (p ms and (1 — (p ms ), 
respectively. 



Equation ( 3.17 ) enables us to obtain a KJMA approx- 
imation for the variance of the Ising magnetization. The 
variance is obtained from the correlation function as 



Var[m(*)] = L^^G^i) . (3.18) 

i=l 

Combining Eqs. (jUjj), (|3T7|), and (EDI), we obtainEl 



G 



L d Var[m(i)] « [m ms - m s ] 2 dU d (2vt) d ^(t) 



9 d (J^ d+1 ) - 9 



iris 



-[1 *bT>£ 



where the function 

e d (x) = 



v 



(3.19) 



(3.20) 



is obtained by numerical integration. Here Xt ^ s the 
isothermal susceptibility in the equilibrium phase, and 
Xt s can be interpreted as an analogous measure of the 
subcritical fluctuations in the metastable phase. 



IV. NUMERICAL RESULTS 

In this Section we present our simulation results for 
the two-dimensional Ising lattice-gas model and use them 
to obtain the parameters in the. theoretical predictions 
of the extended KJMA theoryJjBBEl The relaxation 
func tion, which is a one-point function, is discussed in 
Sec. [VA. Those quantities which also require knowl- 



edge of two-poi nt correlation functions are treated in 
Sees. IVBTVD. In the remainder of this paper we use 
units such that J = ks = 1 and measure time in Monte 
Carlo steps per site (MCSS). 

All the results shown correspond to T = 0.8T C , and 
most of them are for L = 256. Only for the two weakest 
fields, |il| = 0.15 and 0.12, did we use L=1024 to ensure 
a sufficiently large value of L/Rq. Results are averaged 
over 100 independent realizations, except for |if|=0.12, 
for which only 50 realizations were performed due to the 
long lifetime at this weak field. 

The evolution of the system geometry during the de- 
cay process is illustrated in Fig. |l| by a series of typical 
simulation snapshots for |if|=Qrl5. 

Values at 0.8fc B T c ~ 1.815E3 of quantities for the 
Ising model, that are needed to compare the simula- 
tions with the KJMA predictions are as follows. The 
surface tension ctq pi0.745E3 and the zero- field magneti- 
zation m sp «0.954.E£lThe corresponding .jzalues of fid and 
S(T) are fl 2 « 3.153 « tt and 3 « 0.506O This value of 
f^2 gives the constant in the definitions of to and Rq as 
Am 0.870, and it implies that the average deviations of 
the droplets of stable phase from the circular shape as- 
sumed by the KJMA theory should be negligible at this 
temperature. Considering the irregular shapes of the in- 
dividual droplets in Fig. [j], we find it quite remarkable 
that the extended KJMA theory nevertheless gives a very 
good description of the decay process, as we now proceed 
to demonstrate. 



A. Relaxation Function and Metastable 
Magnetization 

Monte Carlo (MC) and fitted KJMA results for the re- 
laxation function <f>(t) are shown together in Fig. || for 
\H\ = 0.2 and 0.4 (both in the MD regime), and 0.8 
(slightly beyond the mean- field spinodal). The results 
are shown on a linear scale in Fig. ||(a) , while the lin ear 
dependence of In (f)(t) on (t/r) 3 predicted by Eq. (3.7) is 
illustrated in Fig. [|(b). The KJMA expression for (f>(t) 
contains two parameters, m ms and Iv 2 , which were de- 
termined by fitting to the M C data as described below. 

As discussed in Sec. IIIB, the theoretical expression for 



(f>(t) follows from the assumption that for times when the 
mean size of the domains of stable phase is much larger 
than the critical droplet size R C (H,T), the mesoscopic 
structure of the Ising model resembles the coa rse-g rained 
KJMA picture. This assumption leads to Eq. ( 3.16| ), from 
which the following two-parameter expression for In <p(t) 
results: 

In [ftp, a(\H\),b(\H\))} a In [a(\H\)} b(\H\) t 3 , (4.1) 

where a(|£f|) contains information about the magnetiza- 
tion of the metastable phase and is given by 



a(\H\ 



m(0) — m s 



and 6(| H |) is obtained from Eq. (3.7) as 



H\H\) 



n 2 iv 2 



(4.2) 



(4.3) 



The theoretical results for <j>(t) were obtained by per- 
forming unweighted linear least-squares fits of Eq. (4.1) 
to the MC results for In <j>{t). Since coalescence effects arc 
expected to make Eq. (4.1) invalid for m < 0, we used 
only data for t <t(\H\) in the fits. We eliminated the 
early-time regime of rapid approach to "metastable equi- 
librium" [the "hooks" most easily seen in Fig. ||(b)] by 
selecting the lower limit of the fitting interval, t m in(|^|)- 
Two different criteria were used. 

(o>) imin(|-ff|) was selected to give a joint extremum for 
a(|_ff|) and b(\H\), yielding lower bounds for m ms and 
Iv 2 . This minimizes the sensitivity of the estimates to 
the cutoff. 

(b) t m i n (\H\) was selected to give a minimum or a plateau 
in the %-square per degree of freedom in the fit. 
For small \H\ the difference between the estimates re- 
sulting from these two criteria is much smaller than their 
individual statistical errors, indicating that the KJMA 
parameters are well defined. For \H\ > 0.3, the time 
scales corresponding to the fast relaxation towards the 
metastable quasi-equilibrium and the slow decay towards 
equilibrium are not well separated. This results in the 
disappearance of the extrema used to select t m in(|-ff|) and 
a rapid loss of precision in the definition of the fitting pa- 
rameters with further increase of \H\. 
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We use unweighted fits because the values of m(t) at 
different t are not independent. The statistical errors in 
the data in the intermediat e-tim e regime expected to be 
most compatible with Eq. (4.1) are larger than those in 
the early- time " hook" regime (see discussion of Var[m(t)] 
in Subsec. IV C[ below). A weighted fitting procedure pro- 
duced much inferior agreement between the fitting func- 
tion and the data than the unweighted procedure. As 
seen from Figs. ||(a) and (b), the agreement between the 
MC results and the predictions of the KJMA theory is 
excellent for intermediate times. 

The progressive breakdown of the validity of the 
KJMA prediction for <fi(t) at late times is illustrated in 
Fig. ||(c), which shows the same data as Fig. ||(b) up 
to t = 2T(\H\). For |ff|=0.2, the KJMA approximation 
agrees very well with the MC data at intermediate times, 
whereas for long times it decays more slowly than the 
MC results. This is expected since the KJMA approxi- 
mation does not incorporate the interface-tension effects 
which accelerate the decay in the late-time regime where 
droplet coalescence becomes important. For \H\ > 0.5, 
the KJMA approximation for <fr(t) agrees well with the 
MC data at intermediate times, whereas for late times 
it decays faster than the MC results, as illustrated by 
the data for \H\ = 0.8. This qualitative change in the 
late-time behavior of 4>(t) signals the breakdown of the 
KJMA nucleation-and-growth picture as \H\ approaches 
the mean-field spinodal (Hmfsp ~ 0.75 at 0.8T C ). In 
the strong-field regime beyond Hmfsp the nucleation of 
very small droplets of stable phase becomes the dominant 
decay mechanism, rather than the growth of larger do- 
mains. The almost perfect agreement between the Ising 
and KJMA results for \H\ =0.4 we believe to be the re- 
sult of an accidental cancellation of corrections at late 
times. 

Monte Carlo data for <fi(t) in the strong-field regime, 
at \H\ = 0.8, 1.0, 2.0, and 3.0, are shown in Fig. [| As 
expected, the MC data are not well approximated by 
the KJMA result in this regime. The solid curve is the 
exact limit for \H\ — > oo; <fr(i) — exp(— t). The data for 
| H | =3.0 are close to this limit. 

For the two weakest fields, \H\ — 0.15 and 0.12, we 
noticed a slight increase in the minimum ^-square per 
degree of freedom in the fits. This may indicate that for 
even weaker fields one may need to consider the "incu- 
bation time" for near-critical cl u&tqrSj which is discussed 
by Shneidman and coworkers .E3EjE3 Investigations for 
weaker fields are t herefore desirable. 

Using Eq. ( |4.2D together with the fitting parameters 
a(|iJ|), we obtained estimates for the metastable mag- 
netization m ms as a function of H. The equilibrium 
magnetizations m s for each value of H, which are nec- 
essary inputs for this calculation, aiere obtained by stan- 
dard equilibrium MC simulationO They are shown in 
the right-hand part of Fig. (^) . 

The estimates for shown in the left-hand part 

of Fig. (||). The statistical errors for these estimates (as 
for most of the other estimates of nonequilibrium quanti- 



ties presented in this paper) were calculated by dividing 
the set of independent runs into five equal batches. They 
are everywhere smaller than the symbol size and there- 
fore not shown. The metastable magnetization is the 
quantity which is most sensitive to the short-time cutoff 
imin(l-HI) used in the fitting process. For \H\ < 0.25, the 
two estimates agree to within the (small) statistical error. 
For stronger fields, the estimates differ, indicating that 
"imsd-ff |) becomes increasingly ill defined as the field is 
increased. This is the main source of uncertainty in our 
estimates of the metastable magnetization. 

The metastable magnetizations shown in Fig. || ap- 
proach the curve of equilibrium magnetizations in a fash- 
ion that appear.s njiiJ/' .smooth and resembles an analytic 
continuationoBrd'Ej We find this resemblance quite re- 
markable, since these estimates are obtained from obser- 
vations of the time-dependence of the decay process, ef- 
fectively using the theoretical KJMA result for <f>(t) to 
subtract those regions of the system which have already 
decayed into the stable phase at any particular time. 
Somewhat fancifully, .this might be called "analytic con- 
tinuation on the fly."E2l 

To check the KJMA estimates for m ms , we choose the 
transfer-matrix (TM) method first suggested and devel- 
oped by Schulman and collaborators.E3~ED A brief de- 
scription of the method with details of its application 
to the present problem is given in Appendix B|. The 
transfer-matrix estimates for m ms , based on Nxoo Ising 
systems with N — 5, . . . , 9, are shown in Fig. [| as solid 
black points. The agreement is gratifying and indicates 
that the metastable order parameter estimates extracted 
from the dynamic MC simulations using the KJMA the- 
ory are consistent with a different theoretical approach 
which is completely independent of the dynamics. 



B. Growth Velocity and Nucleation Rate 

On the basis of the relaxation function alone, one can 
obtain the combination Iv 2 of the nucleation rate and the 
radial growth velocity from the fitting parameter 6(|_ff|). 
To obtain separate estimates for I and v, one needs to 
compare the MC results for the variance of the magneti- 
zation, Var[m(t)], with the corresponding KJMA predic- 
tion given in Eq. ( 3. IE ). If we use the values of m ms and 
Iv 2 obtained from the fits to the relaxation function, and 
Xx as obtained from the fluctuations in the equilibrium 
simulations, then v 2 an d \t s may be determined from 
a linear fit of Eq. ( jn<3| ) to the MC data for Var[m(i)] 
at each value of \H\. For the same reasons discussed 
in the context of the fits to the relaxation function in 



Sec. [V A , we found that an unweighted fitting proce- 
dure was more stable and yielded better overall agree- 
ment than weighted fitting. For each value of |if|, the 
values of Var[m(i)] were fitted over the same time inter- 
val as the corresponding relaxation function, and error 
bars were estimated in the same way as for m ms . An 
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example of such a fit is shown in Fig. |5|. 

In Fig. ^| we show the fitted values of the radial growth 
velocity v(\H\) for \H\ between 0.12 and 0.8. Since this 
method is a rather indirect way to obtain the average ve- 
locity of a convoluted, driven interface, one may reason- 
ably ask whether the fitted is anything more than 
a phenomcnological parameter. In order to answer this 
question, we performed additional numerical and theo- 
retical analyses as discussed in the next two paragraphs. 

In order to further test the identification of v(\H |) with 
an average interface velocity, we performed additional 
MC simulations of the time, evolution of plane interfaces 
driven by an applied field.E2l We started with 64x64 sys- 
tems with all spins antiparallel to the applied field, ex- 
cept for one row of overturned spins along one of the 
lattice edges, and periodic boundary conditions in the 
direction parallel to the resulting interface. We then let 
the systems evolve in an applied fiel d a ccording to the 
Glauber transition probability, Eq. (2.6) with T=0.8T C 
as before, except for the following essential modification. 
Nucleation in the bulk metastable phase was suppressed 
by setting equal to zero the transition probability of any 
spin parallel to all of its nearest neighbors.Ea To distin- 
guish them from the unconstrained growing interfaces 
discussed elsewhere in this paper, we call these inter- 
faces "tame." We performed 100 independent simula- 
tions, continuing each until the interface touched the op- 
posite wall of the simulation box. The average interface 
position in the growth direction was calculated from the 
time-dependent magnetization as y(t) = [m(t) + l]L y /2, 
where L y = 64 is the extent of the simulation lattice in 
the growth direction. Velocity estimates, obtained from 
linear fits to y(t) and averaged over the 100 independent 
runs, are also shown in Fig. ^. They lie very close to a 
straight line through the origin, of slope slightly less than 
that obtained from the fits to the KJMA theory. This is a 
reasonable result, since the absence of subcritical fluctu- 
ations in the "chilled metastable phase" in front of these 
"tame" interfaces should slow down their progress and 
make their average velocity a lower bound for the veloc- 
ity of an interface growing into a metastable phase with 
a thermal distribution of subcritical fluctuations. 

The interface of a growing Ising cluster is in the d#-. 
namic univerpJity class of the Kardar-Parisi-ZhangcJ 
(KPZ) modelEirO The growth velocity of a planar inter- 
face is therefore expected to be linear in H for weak fields, 
as is also predicted— JjorpLarge droplets by the Lifshitz- 
Allen-Cahn theory.EZlL^O However, neither theory ex- 
plicitly gives the proportionality constant, which should 
depend on the average orientation of the interface, as 
well as on the specific, microscopic dynamic. Recently, 
Rikvold and KolesikH have developed an approximate 
theory for the growth velocity of "tame" Ising interfaces, 
based on the solid-on-solid (SQS)-approximation for the 
equilibrium interface structure. This theory gives rise 
to the theoretical curves shown in Fig. |[ The corre- 
sponding analytic expressions are given in Appendix |^. 
The agreement with the MC results for the "tame" in- 



terfaces is excellent. 

The theoretical and numerical evidence presented here 
strongly supports the assertion that v{\H\) obtained by 
fitting KJMA predictions to MC simulations is indeed a 
reasonable estimate of the average propagation velocity 
for the convoluted interfaces which separate the regions of 
stable and metastable phase. The estimate appears satis- 
factory, even though the nucleation rate in the metastable 
phase is too high to measure the growth velocity by more 
direct methods. 

Except at very early tim es, th e magnitude of the term 
proportional to Xt s m Eq- ( {3.1£ ) is much smaller than the 
first term. As a result, the fitted values of v(\H |) are quite 
insensitive to x™ s , which shows large error bars and large 
fluctuations with respect to \H\. Accurate determination 
of Xt s evidently would require much larger data sets than 
used in this study. 

With the separate estimates for b(|-H]) and v(\H\), we 
can easily calculate the nucleation rate I(T, \H\) from 
Eq. (p~g|). The results are shown in Fig. ^. Adjusting the 
unknown coefficient B(T) in the exact asymptotic rela- 
tion Eq. (3.3), such that the theoretical line goes through 
the data point at |i/|=0.15, we find good overall agree- 
ment. The curvature of the theoretical result is due to 
the prefactor exponent K=3. 



C. Correlation Functions 



In Sec. IV B we used Var[m(i)], which is proportional 
to the spatial integral of the correlation function G(fi,t), 
to determine the radial growth velocity v(|i?|). We now 
proceed to obtain the circularly avera ged K JMA corre- 
lation function G(r,t) from Eqs. ( 3.11 ) and ( 3.17 ), using 
the MC estimates for a(\H\), b(\H\), and v(\H\). Since 
the in-phase correlation functions, G ms (r) and G s (r), 
are of very short range, we here set them equal to zero 
for nonzero r. The results for |if|=0.2, 0.4 and 0.8 at 
t—T(\H\) are shown in Fig. ||, together with the corre- 
sponding MC results. The agreement is quite good, ex- 
cept for r w 0. This small discrepancy arises from the 
coarse-grained nature of the KJMA theory and is con- 
sistent with the very short range of the in-phase corre- 
lations. By comparing the theoretical and MC correla- 
tion functions in Fig. pi we infer that the ranges of the 
in-phase correlation functions are on the order of one lat- 
tice constant. The difference between the theoretical and 
simulated correlation functions at r=0 gives an estimate 
of G ms (0) + G s (0). Monte Carlo results for G(r, r(\H\)) 
in the strong-field regime are shown in Fig. |9| for com- 
parison. Note the very short range of the correlations in 
this regime. 

The time evolution of the correlations and the break- 
down of the agreement between the KJMA approxima- 
tion and the MC data for late times and for increasing 
fields are well illustrated by the time-dependent charac- 
teristic length (r(t)), defined in Eqs. (2.11) and (3.13) for 







MC and KJMA, respectively. These quantities are shown 
together in Fig. |o| for (#1=0.2, 0.4 and 0.8. For early 
and intermediate times, the MC and theoretical results 
for (r(t)) increase approximately linearly with t until they 
reach a maximum at a time somewhat beyond r. For late 
times, the MC and theoretical results differ considerably. 
This is easily understood, since the long-time dynamical 
behavior of the Ising model is dominated by interface ten- 
sion effects which are not included in the KJMA theory. 
The characteristic lengths obtained from the MC simu- 
lations are shorter than the KJMA estimates by varying 
amounts, which are less than about 0.5 for t < t(\H\). 
We believe this reflects the short-range, in-phase correla- 
tions. These are ignored in the KJMA estimates, whereas 
they are present in the MC correlation functions, weight- 
ing the latter slightly towards smaller r. 

In Fig. [ll] we show the H dependence of the most 
important lengths that characterize the system during 
the decay. As expected, Rq and the maximum value 
of (r(t)), r max (T, \H\), are proportional over the whole 
range of fields studied. The diameter of a critical droplet, 
2R C (T,H) w a /[m s \H\}, is everywhere smaller than 
these mesoscopic characteristic lengths. 



is not surprising, since it is exactly at these mesoscopic 
length scales that the KJMA theory is expected to de- 
scribe the spatial structure. 

The behavior for large q is best seen in the log-log plots 
in Fig. |l3|(b). The KJMA correlation functions are Ha= 
ear for small r and therefore agree with Porod's Law,E21 
which states that the structure factor for a two-phase 
system with interfaces of negligible thickness should be- 



have as S(q) 



-(d+i) 



for large q. The small oscilla- 



The expression for T(r,t), Eq. (3.11), can be recast in 
the two-parameter scaling form, 



r(r,t)/r(o,t) 



_ f {exp [(r /i? ) 3 * 2 (r/2r )] 



1}, r<2r 
r > 2r Q 

(4.4) 



where ro(t)=vt is proportional the-average radius of the 
growing domains of stable phase. 
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This two-parameter 
scaling behavior is illustrated in Fig. |l^, which shows 
MC and KJMA correlation functions vs r/2ro for two 
different sets of \H\ and t in the MD regime, chosen 
such that they give the same value of tq/Rq. The func- 
tions are normalized such that for both sets of param- 
eters, the KJMA correlation function eq uals unity at 
r=0. The dependence on tq/Rq in Eq. (4.4) may ex- 
plain the breakdown with increasing volume fraction of 
the one-parameter scaling in terms of r/ro, recently used 
by Huang et aln for the experimentally obtained domain 
correlation function of polymer films undergoing phase 
transformation. 



D. Structure Factors 

The good agreement between the MC and KJMA re- 
sults for the correlation functions in the MD regime 
should be accompanied by similar agreement for the 
structure factors. This is confirmed by Fig. 13, which 
shows the MC and KJMA results for the circularly aver- 
aged structure factor S(q,r(\H\)) for |ff|=0.2, 0.4, and 
0.8. 

The extent of the agreement for small q is best seen in 
Fig. Il3](a), which shows the data on a linear scale. This 



tions superimposed on the q~ 3 tails are due to the sharp 
cutoff in r at r=2rg. For the MC data, the thermal 
fluctuations and the lattice cutoff at the Brillouin zone 
boundary causes marked deviations from Porod's law. 
Considerably weaker fields and consequently larger val- 
ues of r max would be necessary to obtain the separation 
of length scales necessary to observe Porod's law in a 
substantial range of q for the MC data. 



V. DISCUSSION 

This paper reports a detailed theoretical and simula- 
tional study of the transient spatial structures that evolve 
during phase transformation driven by a difference in 
free-energy density between a metastable and a stable 
phase. This process is a prototype of metastable decay in 
a wide range of physical and chemical systems, which are 
commonly analyzed in terms of the KJMA theory or one 
of its many extensions and generalizations. The model 
system used in the numerical part of our study is a two- 
dimensional, kinetic Ising lattice-gas model. This is one 
of the first detailed attempts to verify the KJMA theory 
and identify the limitations to its validity, using a model 
system in which the elementary kinetic processes act on 
time and length scales much smaller than those charac- 
teristic of the mesoscopic stable-phase droplets. Since 
the model contains no impurities or free surfaces, the de- 
cay occurs via homogeneous, progressive nucleation and 
subsequent growth of droplets of the stable phase. While 
homogeneous nucleation is less common in nature than 
heterogeneous nucleation at impurities and surfaces, the 
limitations to the KJMA picture that we identify, should 
be valid also for these more complicated situations. 

Our numerical results confirm that the KJMA theory 
for the volume fractions of stable and metastable phase 
(i.e., one-point functions), together with Sekimoto's ex- 
tensions that provide two-point correlation functions, 
give a remarkably accurate description of the decay pro- 
cess for a wide range of system parameters. This regime 
extends surprisingly far towards strong fields (large su- 
persaturations) and the correspondingly small nucleation 
barriers. The conditions for the theory's validity are es- 
sentially as follows. 

• The critical droplets of stable phase must be larger 
than the lattice constant, while at the same time 
much smaller than the system itself. 
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• The system must be sufficiently large that the total 
number of droplets is large. 

• Due to the effects of droplet coalescence, the the- 
ory breaks down for late times, when the remaining 
metastable volume fraction becomes less than one 
half. 

When these conditions are satisfied, there is a large sepa- 
ration between the microscopic time scale and the meso- 
scopic time scale characteristic of the phase separation. 
Under these circumstances, we find the extended KJMA 
theory to be sufficiently accurate that it enables us to 
measure nonequilibrium thermodynamic quantities, in- 
cluding the order parameter in the metastable phase, the 
droplet nucleation rate, and the average propagation ve- 
locity of the convoluted interface between the two phases. 
We are able to verify the measured values by indepen- 
dent theoretical arguments. Due to the relatively short 
lifetime of the metastable phase, these quantities are not 
easy to measure directly by more traditional methods, 
and the methods developed here may therefore be useful 
for experimental measurements on systems undergoing 
phase transformation as well. 

While we have demonstrated excellent agreement be- 
tween the extended KJMA theory and the kinetic Ising 
model in two dimensions and for moderately strong fields, 
we believe it would be useful to perform similar studies 
for weaker fields and. in higher dimensions. A study for 
d = 3 is in progress. t3 
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Si = [m nls + S nls (fi)]u(fi) + [m s + S s (fi)](l - u(n)) , 

(Al) 

where S ms (ri) and 5 s {fi), which both average to zero, 
are local fluctuations in the metastable and the stable 
phase regions, respectively. Assuming the local fluctu- 
ations in the two phases are mutually uncorrelated and 
uncorrelated with the phase field u, Eq. (2.1C) then gives 
Eq. ( |3.17D with 



G ms (r) = (S ms (x)S ms (x + r))(u{x)u{x + r)) (A2) 

and analogously for G s (r). Further assuming that the 
correlation lengths for the local fluctuations are much 
shorter than for the phase field, 

G ms (r) w (S ms (x)6 ma (x + r))(u 2 (x)) 
= (S ms (x)S ms (x + r))(u(x)} 
= (S ms (x)S ma (iS + r))(p ma (A3) 

for all r such that (5 ms (x)S ms (x + r)) is nonzero. As a 
result, the spatial sum over G ms (r) can be considered 
an "isothermal susceptibility" for the metastable phase, 
weighted by the metastable volume fraction: 



]TG ms (r 



(A4) 



i=i 



The same reasoning is applied t o the fluctuations in the 
equilibrium phase, yielding Eq. (3.19). 
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APPENDIX A: COARSE-GRAINED 
APPROXIMATION FOR THE CORRELATION 
FUNCTIONS 

The coarse-grai ned a pprox imat ions for G(f,t) and 
L d Va,r[m(t)}, Eqs. ([3.17D and (|3.19| ), are obtained as fol- 
lows (with the time variable suppressed for simplicity of 
notation). We write the local spin variables as 



APPENDIX B: TRANSFER-MATRIX 
ESTIMATES OF THE METASTABLE 
MAGNETIZATION 

Here we summarize the transfer-matrix (TM) calcula- 
tion used to check t he con sistency of the KJMA estimates 
for m ms (H) in Sec. IV A. 



The field-theoreti cal r esult for the nucleation rate used 
in this work, Eq. ( |3.3| ), is proportional to the imagi- 
nary part of an analytic continuation, qf t he equilibrium 
free energy intoihc metastable phase-ESO Schulman and 
collaboratorsL3lE3 suggested that this analytic continua- 
tion could be found from some of the subdominant eigen- 
values and eigenvectors of the TM commonly used to 
calculate the equilibrium partition function rasd correla- 
tion lengths of Ising and lattice-gas systems,t22l and they 



successfully obtained the metastable .order parameter for 
the Ising ferromagnet at T w 0.4T c .l_a The method has 
later been extend ed to obtain the full, complex-valued 
free energy.oEaO 

The dominant eigenvalue Ao of the TM for an N x 
oo Ising system is positive—and nondegenerate by the 
Perron- Frobenius theorem. E£5l It is related to the free 
energy per site as F(T,H) — ~(T/N)\nX , and the 
equilibrium magnetization is obtained as m eq (T, H) = 
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—dF(T, H)/dH = (0\M\0), where M is the magnetiza- 
tion operator and the bra and ket are the left and right 
eigenvectors corresponding to Ao- When the subdomi- 
nant eigenvalues A a are plotted in the same logarithmic 
way used to obtain the equilibrium free energy from Ao, 
one obtains a plot like the one shown in Fig. With 
each branch in the figure one can associate a magne- 
tization m a (T,H) = d[(T/N)bx\X a \]/dH = (a\M\a). 
The metastable branch corresponds to the union of the 
lowest-lying eigenvalue branches in the figure, which have 
a magnetization whose sign is opposite that of the equi- 
librium magnetization. It is marked with thick curve 
segments in Fig. [l4|. At specific values of H the eigenval- 
ues in the composite metastable branch undergo avoided 
crossings with other branches, at which their eigenvectors 
and magnetizations vary rapidly. The field correspond- 
ing to the nth crossing (n = corresponds to Hq = 0) 
depends on n and the strip width N approximately as 
H n w 2M(N — n), with best agreement for low T and 
small n.c3 The magnetization along the branch between 
the (n— l)th and nth crossing is what we refer to as the 
"nth lobe" in the caption of Fig. [|. The point on each 
lobe, where the magnetization calculated from that lobe 
has its extremum, is marked as a solid black circle in 
Fig. [l4| These points correspond to the points similarly 
marked in Figs. |J and [H]. 

For this work we numerically diagonalizcd-thc TM for 
N = 5, . . . , 9, using the subroutine jacobiLH in double 
precision on a DEC- alpha workstation and a Cray J90 
supercomputer. As seen in Fig. |J, we found excellent 
agreement between the KJMA results and the TM esti- 
mates from the first lobe for very weak fields and from 
the second lobe for stronger fields. The values of m ms (7J) 
extracted from the second lobe for the different values of 
N agree with the KJMA estimates to within the uncer- 
tainty in the latter. Estimates based on the third and 
higher lobes give less satisfactory agreement. 



APPENDIX C: APPROXIMATE EXPRESSION 
FOR THE INTERFACE VELOCITY 

The analytic approximation used here for the propaga- 
tion velocity of a field-driven "tame" interface in a two- 
dimensional, square-lattice Ising model with-Glauber dy- 
namics will be described in detail elsewhereJ£3 It uses the 
SOS approximation for the equilibrium interface struc- 
ture (which remarkably gives the exact surface tension 
for interfaces parallel to the lattice directions-and an ex- 
cellent approximation for inclined interfacesL^IEj) to esti- 
mate the populations in the different spin classes used-La 
the n-fold way rejection-free Monte Carlo algorithm.EH 
These class populations are used together with the con- 
tributions to the average propagation velocity from spins 
in each class, which are easily obtained from the transi- 
tion rates, to calculate the overall average velocity. For 
the special case of Glauber transition rates, isotropic in- 
teractions, and an interface which is on average parallel 



to-one of the symmetry directions of the lattice, the result 



(.(T.iO>-^{« 



1 + X 2 



sinh(2/3J) 
cosh(/3H) 



X 2 



1 - A 2 



X 2 



1 



2(1 + 2X) 

( S inh(2/3 J) \ ' 
' \cosh(l3H) J 



(CI) 



Here X = exp(— 2 (3 J) corresponds to a linear-response 
like approximation, in which the average class popula- 
tions are approximated by their equilibrium values at 
H = 0, whereas X = exp(— 2/3 J) cosh(/3H) yields a 
nonlinear-response approximation which accounts for ef- 
fects of the applied field on the nonequilibrium class pop- 
ulations. In Fig. ||, (v(T,H)) from Eq. (CI) is shown 
vs | if | at T = 0.8T C . The agreement of the nonlinear- 
response result with the directly simulated "tame" inter- 
face velocities is excellent. Since the surface tension at 
0.8T C is very close to isotropic, results for inclined inter- 
faces are not needed here. 
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FIG. 1. Snapshots showing the time evolution of the sys- 
tem configuration during a typical Monte Carlo (MC) simula- 
tion for L=250, T=0.8T C , and |ii"|=0.15, for which the lifetime 
r ~ 392 MCSS. Light gray represents the metastable spin di- 
rection and black the stable spin direction. While supercriti- 
cal droplets in one snapshot can also be identified in the next, 
the subcritical fluctuations are uncorrelated between snap- 
shots. The two characteristic lengths, Ro ~ 25 and R c ~ 2.5, 
are shown by the long and short bar below each snapshot, 
respectively, (a) t = 80 MCSS. A few supercritical droplets 
are seen, but most of the black pixels represent subcritical 
fluctuations, (b) t — 260 MCSS. Many large supercritical 
droplets are growing against the metastable background, (c) 
t = 390 MCSSs; t. The magnetization is close to zero, sev- 
eral of the supercritical droplets have coalesced, and the stable 
phase is close to percolating. Microscopic equilibrium fluctu- 
ations appear as light specks inside the stable-phase regions. 
Figure courtesy of G. Korniss. 
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FIG. 2. MC (points) and fitted KJMA results (lines) for 
the relaxation function of the Ising model, 4>(t), for L=256 
and T=0.8T C . Results are shown for |#|=0.2 (O and solid 
curve) and 0.4 (□ and dot-dashed curve), both in the MD 
regime, and for |J/|=0.8 (O and dashed curve), slightly be- 
yond the mean-field spinodal. The lifetimes are r ~ 186 
MCSS for |#|=0.2, r » 41 MCSS for |#|=0.4, and r » 12 
MCSS for |_H~|=0.8. (a) Linear scale vs t. Derivatives of 
data such as these pro-wide estimates for transient currents 
in ferroelectric. pyitchingtJ and electrochemical potential-step 
experiments.LJ'Cj (b) Semi-logarithmic scale vs (t/r) 3 . For 
clarity, the data for |_ff|=0.8 and 0.2 have been displaced by 
±0.2, respectively, (c) Same as (b), but plotted for t/r < 2 to 
show the expected deviations of the MC data from the KJMA 
approximation for later times. 
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FIG. 3. MC results for (j>(t) in the strong-field regime 
at |-ff|=0.8 (o and dashed curve, repeated from Fig. H for 
comparison), 1.0 (O), 2.0 (□), and 3.0 (A). The lifetimes are 
t « 8.42 MCSS for [H|=1.0, r » 2.70 MCSS for |ff|=2.0, 
and r fa 1.42 MCSS for \H\=3.0. The solid curves represent 
exponential relaxation, which is the exact result in the limit 
\H\ —* oo. (a) Linear scale, (b) \n<f>(t) vs t, emphasizing the 
approach towards exponential decay with increasing \H\. 
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FIG. 4. Stable magnetization, m s (O), and metastable 
magnetization, m ms (o and □), shown vs H. The data points 
for m ms re presen ted as o are based on the selection criterion 
(a) in Sec. 1VA, and those represented as □ on criterion (b). 
The data are for L = 256 except for H = —0.15 and —0.12, 
which correspond to L = 1024. The thin solid curves are trans- 
fer-matrix (TM) results. For H > the curve represents the 
equilibrium magnetization as obtained for an N x oo system 
with N — 9, which is seen to be in complete agreement with 
the MC results. For —0.2 < H < the curve represents the 
metastable magnetization corresponding to "the first lobe" 
[see Appendix j3[ for iV = 9, while the solid circles represent 
the maximum of this lobe for N = 5, . . . , 9 from left to right. 
For H < — 0.2 the curves and solid circles represent "the sec- 
ond lobe" for N — 5, . . . , 9 from left to right. 
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FIG. 5. An example of L 2 Var[m(£)] as obtained from 
a comparison of 100 independent MC simulations (O) and 
from a lea st-squares fit of the KJMA theoretical expression, 
Eq. (3.19), (solid curve). This particular result corresponds 
to | jy |=0.2. The thin vertical lines mark the time interval 
[imm,i~], over which the fitting was performed. 




FIG. 7. The nucleation rate I(T, \H\), computed using 
Eq. (4.3) with values of b(\H\) and obtained from fits 

to 4>(t) and Var[m(i)], respectively. The dashed curve is a 



one-parameter fit to the exact asymptotic result, Eq. (3.3) 
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FIG. 6. The average radial growth velocity of the domains 
of stable phase, v(T,\H\), vs \H\. The circles and squares 
(corresponding to the same selection criter ia as in Fig. ft) are 
the velocities obtained by fitting Eq. (3.19) to the MC results 
for Var[m(t)]. The diamonds are the velocities obtained in 
MC simulations of a growing plane interface in which nucle- 
ation inside the single-phase domains is suppressed - a "tame" 
interface. The thin dashed curve and the solid curve represent 
the analytical linear-response and nonlinear -i^espon se ap prox- 
imations for the "tame" interface velocitiesJ23 Eq. (CI). 
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FIG. 9. MC results for the correlation function 

G{r,r(\H\)) in the strong-field regime for \H\=1.0 (O), 
|Jf|=2.0 (□), and |ff|=3.0 (A). The lines connecting the 
data points are merely guides to the eye. Note the very short 
ranges of these correlation functions. 
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FIG. 10. MC and KJMA results for the first moment of 
G(r,t), (r(t)). Note the reasonable agreement at early and 
intermediate times. The rapid decrease of the MC character- 
istic length at late times reflects the acceleration of the decay 
of the metastable phase due to interface-tension effects during 
the droplet-coalescence regime. In each panel the thin vertical 
line marks t(\H\). (a) |#|=0.2. (b) |H|=0.4. (c) |#j=0.8. 
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FIG. 11. Field dependence of the characteristic lengths: 
the critical droplet diameter 2R C (solid curve) , the mean sepa- 
ration between droplets of stable phase Ro (O, A, and dashed 
curve) and twice the maximum of (r(t)), 2r max (o). Two sets 
of estimates for R$ are shown: Av/(iw 2 ) 1/ ' 3 (O with error 
bars) and vr (A and dashed curve). The lengths that charac- 
terize the mesoscopic structure, _Ro and r max , remain propor- 
tional over the whole range of fields studied. The two chains of 
solid circles relate the strip width s N of the transfer-matrices 
used to calculate m ms (Sec. IV A and Appendix |§|) to the val- 
ues of H for which m ma so calculated has a maximum. The 
data points correspond to those shown by the same symbols 
in Fig. The chain between H = and 0.1 represents the 
first lobe, and the chain between H = 0.2 and 0.6 represents 
the second lobe; in both cases N — 9, 5 from left to right. 
Comparison of N to the characteristic lengths indicate that 
the first lobe samples only subcritical fluctuations, while the 
second lobe also samples rare supercritical fluctuations in the 
constrained ensemble represented by the TM eigenspaces. 
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FIG. 12. Two-parameter scaling plot of the normalized 
correlation function G (r, r(\H\)) /Gkjma (0, t(\H\)). The 
data points are MC results and the solid curve is the scal- 
ing form of the KJMA correlation function T(r,t), Eq. (4.4). 
The MC results are shown for \H\=0.2 and t=185 MCSS (O), 
and for |#|=0.25 and t=110 MCSS (□). These fields and 
times are chosen such that ro/i?o~0.97 in both cases. 




FIG. 13. The structure factor S (q,r(\H\)). The data 
points correspond to MC results and the curves to theoreti- 
cal results calculated by taking the Fourier transform of the 
KJMA results for G(r,t). Results are shown for |ff|=0.2 (O, 
solid curve), |H|=0.4 (□, dotdashed curve), and |fl"|=0.8 (o, 
dashed curve), (a) Linear scale, (b) Log-log scale. For clar- 
ity, the results for \H\—0.2 and 0.8 have been offset by ± two 
decades, respectively. 
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FIG. 14. The nineteen largest transfer-matrix eigenvalues 
A aj plotted vs H as —(T/N) In |A Q for N = 6. The low- 
est-lying branch (dot-dashed line), which corresponds to the 
dominant eigenvalue Ao, is the equilibrium free energy per 
spin. The metastable branch is represented by the heavy 
curve segments. The solid circles represent the points along 
each segment, where the magnetization m a (T,H) has its ex- 
tremum. These points correspond to those similarly marked 
in Figs. ^ and y. 
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